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Abstract 

We propose an efficient stochastic method to implement numerically the Bogolubov 
approach to study finite-temperature Bose-Einstein condensates. Our method is based on 
the Wigner representation of the density matrix describing the non condensed modes and 
a Brownian motion simulation to sample the Wigner distribution at thermal equilibrium. 
Allowing to sample any density operator Gaussian in the field variables, our method is 
very general and it applies both to the Bogolubov and to the Hartree-Fock Bogolubov 
approach, in the equilibrium case as well as in the time-dependent case. We think that our 
method can be useful to study trapped Bose-Einstein condensates in two or three spatial 
dimensions without rotational symmetry properties, as in the case of condensates with 
vortices, where the traditional Bogolubov approach is difficult to implement numerically 
due to the need to diagonalize very big matrices. 

1 Introduction 

Since the first experimental demonstrations of Bose-Einstein condensation in atomic gases 
^ 0, 1^, we have assisted to a fantastic development of the field both on experimental and 
theoretical ground. From the theoretical point of view the Gross-Pitaevskii equation (GPE) 
or non-linear Schrodinger equation for the condensate wave function, in its steady-state and 
time dependent versions, revealed itself to be an extremely powerful tool to explain most of the 
experimental results p. 

Generally speaking the Gross-Pitaevskii equation is obtained by neglecting the interaction 
between the condensate and non condensed particles [|]. This approximation is justified at 
temperatures much below the transition temperature where the non condensed cloud is very 
small and much less dense that the condensate. This is often the case in experimental conditions 
where almost pure condensates can be obtained. 

There are nevertheless situations in which the rich physics coming from the non condensed 
fraction of the cloud and its coupling to the condensate cannot be neglected. A first example 
concerns the collective modes of the gas induced by a modulation of the trapping potential 
10, Different theories going beyond the Gross-Pitaevskii equation were put forward in order 



to explain the temperature dependence of the frequencies and the damping rates measured in 
experiments ||, |1^, [Tl| (for more references see 0, pages 489-500). 

A second example where the non condensed fraction may play a significant role concerns 
condensates in thermodynamically unstable states, like vortices and dark solitons recently cre- 
ated experimentally [12, |14[ 15]. In |13] for example one or more vortices are created by 



stirring the condensate with a slightly anisotropic rotating potential: vortices appear provided 
that the rotation frequency exceeds a critical value. If the rotation is then stopped the vortices 
are not any more a local minimum of energy and are expected to spiral out of the conden- 
sate [l^, 0]. It is indeed found experimentally that the vortices leave the condensate [|l^]. It 



is predicted theoretically that the dynamics and lifetime of the thermodynamically unstable 
Bose-Einstein condensate depend crucially on the interaction between the condensate and the 
non condensed cloud UlSj . 

The examples we have given above concern time dependent properties of the condensate. 
Another class of interesting finite temperature issues is given by the thermal equilibrium prop- 
erties of the gas. To cite some example, one might be interested in the non condensed spatial 
density of the atoms, or in the thermal fluctuations of the number of particles in the conden- 
sate [T^, Condensates with vortices give rise again to interesting questions in this 
respect such as the thermal fluctuations in the position of the vortex core and the temperature 
dependence of the critical rotation velocity ||2^ . 

A well established technique to study Bose-Einstein condensates at thermal equilibrium is 
the Bogolubov method, which is a perturbative method valid for low temperatures when the 
density of condensate particles is much larger than the density of non condensed particles. This 
method relies on the quadratization of the Hamiltonian in the field operator representing the 
non condensed fraction, and on the Bogolubov transform |]2^. Another method which could be 
valid for a larger range of temperatures is the Hartree-Fock-Bogolubov (HFB) method p5|, ^ , 
which is instead a variational method, based on a Gaussian Ansatz for the system density 
operator. For time dependent problems, one disposes of the former two methods, plus more 
refined techniques including kinetic equations for the non condensed particles [^, |2^, |2^, |5D|. 

A common point to all these methods is the necessity to calculate mean values using a density 
operator which is an exponential of a quadratic form in the field operator: the exponential of 
the quadratized Hamiltonian in the case of the Bogolubov approach or the Gaussian Ansatz for 
the HFB approach. The procedure that is usually followed to perform those calculations is to 
diagonalize the quadratic form by the Bogolubov transform, turning it into a sum of decoupled 
harmonic oscillator energy terms. 

Unfortunately the corresponding calculations are quite heavy for 3D non-homogeneous sys- 
tems such as trapped gases, and this constitutes a serious limitation to the use of these methods. 
E.g. in practice, by discretizing the real space on a grid with a modest number of 64 points per 
spatial dimension, and in the absence of rotational symmetry (as in a multi-vortex configura- 
tion) one has to diagonalize a matrix which is [2(64)" x 2(64)"'] large, where n = 1,2,3 is the 
dimension of the space. For n = 3 the matrix to diagonalize is 524 288 x 524 288. 

In this paper we propose an alternative formulation of the Bogolubov theory which can be 
implemented numerically more efficiently in two or three spatial dimensions since it avoids the 
diagonalization of big matrices. Our method relies on the use of the Wigner representation 
of the density operator and on a Monte Carlo simulation to sample the steady-state Wigner 
distribution of the system at finite temperature. Since our method allows to sample a "Gaus- 
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sian" density operator, it applies both to the Bogolubov and to the HFB theory, in the thermal 
equilibrium as well as in the time dependent case. We give in this paper a detailed derivation 
of the method in the frame of a particle number conserving Bogolubov theory [^, |3^, Q . We 



consider first the thermal equilibrium case and then we show how to extend the treatment to 
the time dependent case. 

The paper is organized as follows: in section ^ we recall the number-conserving version of 
the Bogolubov theory put forward in which we adopt in the following. In section ^ we 



reformulate the theory by expressing the Bogolubov density operator in terms of the Wigner 
quasi-distribution function. In section H we describe the stochastic simulation that we use to 
sample the equilibrium Wigner function. In section |^ we compare numerically our stochastic 
method to the direct diagonalization used in the traditional Bogolubov approach in the case 
of a ID trapped Bose gas, where the explicit diagonalization is feasible, and we calculate the 
spatial distribution of the non condensed fraction. In section ^ we show how a time dependent 
problem would be treated with our method. Conclusions and perspectives are presented in 
section |^. 

2 Bogolubov theory conserving the total number of par- 
ticles 

In this section we recall briefly the main ideas and results of where a number-conserving 
version of the Bogolubov theory is established. The reader already familiar with this theory 
can jump directly to section |. 

We consider a cloud of atoms trapped in the potential U (f) interacting through an effective 
low energy contact potential V = gS{fl — ): 

H = j d^f^\f )hitp{f) + ^ij^{f)ij^{f)ij{f)ij{r) (1) 
where ip is the atomic field operator, hi is the one-body Hamiltonian: 

= ~A + U(r) , (2) 

The coupling constant g is expressed in terms of the s-wave scattering length a of the true 
interaction potential: 

9 = • (3) 

m 

Let pi be the one-body density operator of our system and let 0ex be the eigenvector of pi 
with the largest eigenvalue Nq: 

Pll^e.) = Nolcpe.) . (4) 

We normalize to unity. Physically is the most populated single particle state. We assume 
here the presence of a Bose-Einstein condensate so that (j)cx is the condensate wave function 
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and No ~ A^. This motivates the sphtting of the atomic field operator ip into condensate and 
non condensed modes: 

i^{r) = 0ex(r )a<^,, + 6ip{f) (5) 

where a^^^ annihilates a particle in the condensate wavefunction 0ex; and Siplf) is orthogonal 
to (pe^: 

' €M = . (6) 

An important property that comes from (^) and it that there is no single particle coherence 
between the condensate and the non condensed modes: 

{alM = . (7) 

The interest of the splitting is that the typical matrix element of d^^^ is of order Nq^"^ ^ N^^"^ 
while the matrix elements of 5%jj scale as (SNy^'^ <^ N^^"^ where 5N is the number of non 
condensed particles. 

An essential ingredient of the number conserving Bogolubov approach is the operator Aex 
transferring one non condensed particle into the condensate: 

Aex = N-'/' 6i> . (8) 

This operator plays here the role played by 6iJj in the symmetry breaking approach commonly 
adopted. In particular due to the expectation value of Aex vanishes. 
In the weakly interacting Bose gas limit 

^ oo and 9 N = const , (9) 

and for a constant trapping potential, the mean interaction energy per particle ng (where n 
is the density) tends to a finite value whereas the small gaseous parameter (na^)^/^ tends to 
zero. If the temperature, the trapping potential and the mean interaction energy are fixed, the 
number of non condensed particles 6N is bounded from above so that SN ^ in the large 
limit H. One can then make a systematic expansion of the exact condensate wavefunction 0ex 
and of the fields Aex and AJ^ in powers of the small parameter: 

e = ^6N/N ~ 1/Vn. (10) 

Formally: 

Aex = A(o) + -^A(i) + i^A(2) + . . . (11) 



(0) + ^'^(i) + ^<^(2) + • ■ ■ (12) 



One then calculates d{Acx)/dt = order by order in e. At the order —1 in e (the lowest 
approximation order), one finds A(o) = and 0(o) = (f) where (p is solution of the time- dependent 
Gross-Pitaevskii equation: 

tndt(j)=[h + Ng\<j)\^-fi]<f). (13) 
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The quantity fi in ( ]TB| ) is such that the wavefunction of the condensate at equihbrium is a time 
independent solution of the above equation. Physically fi is the lowest order approximation to 
the chemical potential of the gas. 

The order in e provides the equations of motion for A = --^A(i) and A''': 

*«'(|.)-^(|.)- 

The first non zero correction to 0(o) is j^<P{2)i whose equation of motion is obtained in the next 
order of the expansion. The operator C in ([T^) is given by: 



H^^-^Ji + QNg\(t)\''Q QNg<f)^Q* 
C=\ I (15) 

-Q*Ng{<p*)'Q - [H;^ -fi + Q*Ng\<p\'Q*] 

where 

H^ = h + Ng\(P\^ (16) 

and Q projects orthogonally to 0: 

g = i-|0)(0|. (17) 

The complex conjugation indicated by the star is taken in the r-basis so that e.g. Q* projects 
orthogonally to the state with wavefunction 0*. The projection operator appears as well in the 
commutation relations of A and A"'": 

[A{r),^Ks)] = {r\Q\s) . (18) 

The linear equations of motion ([1^ for A correspond to an approximation of the Hamil- 
tonian which is quadratic in A. This quadratic approximation can be obtained by inserting 
the decomposition (^ in the Hamiltonian, by neglecting terms cubic or quartic in and by 
using: 

«L«<^ex = N-5N (19) 
6N = [ 5i)^5i) ~ /" AU . (20) 



As equation (O) is satisfied, no term linear in A is left in the Hamiltonian. One is left with 



the quadratic Hamiltonian [22 



A 

At 



(21) 



At thermal equilibrium in the canonical ensemble, the non condensed atoms are then de- 
scribed by the density operator: 

^nc = ie-^^-^<^ . (22) 
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The usual way of proceeding at this point is to diagonahze the operator L. Let us introduce 
the eigenvalues and eigenvectors of £: 



We suppose that all the eigenvalues are real so that the equilibrium condensate wavefunction 
is a dynamically stable solution of the Gross-Pitaevskii equation. Moreover we restrict the 
notation {uk,Vk) to the modes satisfying the normalization and orthogonality conditions: 

^fcWfc' - ^fc^fc' = ^k,k' ■ (24) 

To each mode {uk,Vk) is associated a "dual mode" {vl,ul), of energy — and satisfying obvi- 
ously / vlvk' — u\uk' = —Sk^k'- In general the (Mfc,ffc)'s plus the dual modes form a complete 
family in the subspace orthogonal to (0, 0) and (0, (p*) so that the field A can then be expressed 
as: 

^ )-^bk( ) "^^^^ ) ^ (25) 

The normalization condition of the {uk, ffc)'s ensures that the bk, b^ satisfy bosonic commutation 
relations. If (|25| ) is injected in (|2TD, the quadratic Hamiltonian takes the form: 

H^u.d = Eo + Y,^kblbk (26) 

k 

corresponding to a set of decoupled harmonic oscillators. For thermodynamical stability of the 
condensate we must require that all the are positive. Physical quantities can then be readily 
extracted from a^c- The inconvenient of this method is that the diagonalization of C in two 
and three dimensions is a very heavy numerical task in the absence of particular symmetries of 
the problem. 



3 Formulation of the Bogolubov theory in terms of the 
Wigner function 

We wish to deal with complex numbers and complex functions instead of operators and fields 



?fc , Ok 



h* 
Ok 



A\r) , A(f) 



A*(f) , A(f ) . 



For this reason we introduce the Wigner quasi distribution function W{{bk}, {bl}) defined as 
the Fourier transform of the characteristic function x({afc}): 



X({«fc},{afc}) = Tr 



d-nc Y\_ exp {akbl - albk) 



(27) 
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and 



W{{h}, m) = I xi{a,h {al}) Hexp [{alb, - a,bl)] (^^^ "^K^^"^ . (28) 

With the diagonahzed form of the quadratic Hamiltonian (p6D, the density operator ( p2D cor- 
responds to decoupled harmonic oscillators at thermal equilibrium, so that: 



Wi{bk},{bl}) = n2tanh (^^^ exp [-2|6,|2tanh(/5efc/2)] . 



(29) 



In the high temperature limit ksT ^ the Wigner distribution for mode k is simply the 
Boltzmann distribution for a classical field. In the low temperature limit fc^T -C it is worth 
noting that the Wigner distribution has a finite width tending to 1/2: even at zero temperature, 
the c-number quantities bk fluctuate mimicking the quantum fluctuations of the operators bk- 
Since we want to avoid the diagonalization of the operator C, we wish now to express the 
Wigner function (PD| ) in terms of the functions A and A* corresponding to the field operators 
A and A"^: 



A*(f) J ~ ^ V Vk{r) J uUr 

k 

By using ( pOf ) and the properties (0) of the Uk{r) and Vk{r), one can show that the Wigner 
function takes the form: 



iy(A,A*) = det ■ 



Ari tanh 

' 2 



exp|-|(A*,A)- r/tanh^ ( A*) } ^^^^ 



where t] is the matrix: 



"H -Id !• 



and the determinant is defined in the space orthogonal to (0, 0) and (0, 0*). 

The quantum mechanical mean values of totally symmetrized operators are then simply 
given as averages over the Wigner distribution p4|; e.g.: 

(A(rl)A*(r-2))^ = i(At(rl )A(r-2 ) + A(rl)At(rl )). (33) 



Note that according to ([TsD the previous expression gives infinity for rl = r^! In practice we 
will not encounter this problem because we will work on a grid with a finite number of modes. 



4 Brownian motion simulation 

Since we quadratized the Hamiltonian with respect to A and A"!", the Wigner distribution ( PT| ) 
is Gaussian and can be recovered as the steady-state probability distribution of a conveniently 
chosen Brownian motion of the functions A and A*. In this section we show that one can find 
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such Brownian motion in a convenient way for numerical simulations. First of all we discretize 
our problem on a grid. To this aim we introduce the vector A representing the set of values 
{A(r/ )} where = (xj, yj, Zk) are points of the grid: 

(A). = A(r-I). (34) 

If dV is the volume element of the grid, the Wigner function (^Tj) takes the form: 

/ A 



W{k,K*) =^exp 



-dV{K\k)- T^tanh :^ 



(35) 



where now C is represented by a matrix, and where we did not write explicitly the normalization 
factor. The vector A is then submitted to the following stochastic evolution during the time 
step dt: 

dA = Fdt + Yn d^ + Yi2 dC (36) 

where F is a linear friction force 

F = -an A - A* (37) 

and d^ is a noise term satisfying: 



dm = (38) 
dt 

dan)de{r2) = 2— [5.-.„^,-rfV0(fi)0*(f2)] (39) 

dm)dm) = 0, (40) 



where the overline denotes the average over the time interval dt. In (^6]) and (|37D, A and ^ 
are vectors whose length A/" is the size of the grid, while aij and Yij = 1, 2) are JV x JV 
matrices. The prescription (pOf ) is equivalent to the usual prescription for spatially 5-correlated 
Gaussian white noise with the additional requirement that the noise should be orthogonal to 
the condensate wavefunction 0, so that A remains within the subspace orthogonal to 0. It is 
important to note that the choice of the 'time' variable t in the Brownian motion simulation is 
totally arbitrary; in what follows dt will have the dimension of an energy. 

Let P{A,A*,t) be the probability distribution of A and A* at time t. The stochastic equa- 
tions ( PBD can then be turned into a Fokker-Plank equation for P(A, A*,t): 

dtP + [dAif)Ar) + c.c] = (41) 

f 

where J{r) is a density current: 

J(f ) = F(f )P -^Yl dAis)iD,2ir, s)P) + 9A.(,-.)(Z}n(r, s)P) (42) 
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and Dij are constant diffusion matrices: 



dV 



{DnW,s) = —dA{f)dA*{s) 



dV 
2dt 



dA{f)dA{s" 



(43) 
(44) 



Neglecting terms on the order of dt, we can express tlie identities ( ^3| - PD more syntlietically in 
terms of a diffusion matrix: 



D 



wfiere we liave introduced tlie noise matrix 



Y 



^11 ^12 
V* V* 



(45) 



(46) 



Tlie steady-state solution of (^if) with J(r ) = is found to be: 

P oc exp 

with the friction matrix a defined as: 



J^{A\A)-D-^a 



A 
A* 



(47) 



a 



an "12 



a 



12 



a 



11 



(48) 



If we now require that the steady-state distribution ( ^Tj) coincides with the Wigner distribution 
([3T| ) describing thermal equilibrium, we obtain: 



a = 2 Drj tanh 



(49) 



Equation (^) is simply the generalization to a field of Einstein's relation for Brownian motion 
of a particle: in the high temperature limit the diffusion term divided by the friction term is 
proportional to temperature. 

Provided that the Einstein's relation (^9]) is satisfied, the choice of a and Y (which deter- 
mines D) is not unique. It turns out that a convenient choice of the friction and noise terms 
from the point of view of a numerical simulation is given by: 



a 



cosh ( ^— 
2 



sinh ( ^— 
2 



/if 



and: 



Y 



cosh ( 



(50) 



(51) 
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One can check that the choice (^) and ( pT| ) satisfies (^) by using the following properties of 
the operators L and ry: 

T] Lri = C) and 77^^ = 77. (52) 

One also has to check that the choice ( ^0] ) and (^) reproduces the structure of the matrices 



(|46| ) and (0), that is the second line of the matrix is obtained from the first line by complex 



conjugation and exchange of indices 1 and 2, a property that can be formalized as 

SaS = a* and SYS = Y (53) 

where S is the matrix 

Id o)- («^) 

From the properties S CS = —C* and S r] S = —77 one finds that (^) and ( |FI| ) have the right 
structure indeed. 

The convenience of the choice ( |50D and (0) relies in the fact that one then is brought to 
calculate the action of the operators exp on the functions A, which is easily implemented by 
using a splitting plus Fourier transform technique (more details about the numerical simulation 
will be given in the end of the next section). 



5 Test of the method for a ID trapped Bose gas 



As a test of the method we have made a numerical comparison between the Monte Carlo 
simulation and the direct diagonalization of the operator C for the case of a Bose-Einstein 
condensate in a harmonic potential with frequency u in one spatial dimension. Physically 
this situation would correspond to a very elongated cigar shaped trap with a longitudinal trap 
frequency u and a strong transversal confinement huj± > ksT, so that most of the atoms are 
in the transverse ground state of the trap. In these conditions the field operator can be 
approximated by: 



^(f) = (j)^{x,y)i){z) 



(55) 



where the operators 'ip'^{z) and 'ip{z) have the usual bosonic commutation relations for field 
operators in one dimension. We have then reduced our problem to one spatial dimension with 
an effective coupling constant qid [03: 



giD= g dx dy\(j)±{x,y)\ 



(56) 



We wish to calculate the spatial density of non condensed atoms. In order to calculate the 
average {A'^ {z)A{z)) we use the formula (^) and equation (|18|) in its discretized version, to 
obtain: 



(At(^)A(^)) = {A*iz)A{z)) - 



1 



1 

dz 



(57) 
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The result for the non condensed cloud spatial density as a function of z is reported in figure ^ 
We have chosen a total number of atoms N = 10^, a temperature ksT = SOhuj, and a coupling 
constant gm = O.Olh(huj/mY^'^ leading to a chemical potential /i = lA.lhto. The full line is the 
direct diagonalization result and the points with the error bars are obtained with our Monte 
Carlo simulation with 2000 realizations. The non condensed cloud is repelled by the condensate 
which is sitting in the center of the harmonic trap, which explains the double peak structure 
of figure |I[ 

We can also get information on the number of condensate particles by calculating the number 
of non-condensed particles, since the total number of particles is fixed. The mean number of 
non-condensed particles is given by the spatial integral of Eg . (^7|) : 



{6N) = {5N)w - (58) 



where the operator 6N, defined in (pO]), is replaced at the present order by its approximation 
in terms of A, A/" is the number of modes in the simulation (that is the number of points on 
the spatial grid) and 

6N = dzJ2^*i^)Hz)- (59) 

z 

After some algebra, involving the total symmetrization of a product of 4 operators A or A''', we 
also obtain 

{SN') - {6N)' = {6N')w - ml - (60) 
For the parameters of the figure we get the mean values and the standard deviations 

((5iV))Mc = 385±6 and ((5iV))diag = 391 , (61) 

a{5N)Mc = 269 ± 10 and a((5iV)diag = 279 , (62) 

(63) 

where the subscript 'MC represents the results of the Monte Carlo simulation and 'diag' indi- 
cates the predictions of a direct diagonalization of C 

We give in the following some details and tricks about the numerical simulation. The 
first thing do to is to solve the steady-state Gross-Pitaevskii equation to find the function 0. 
This is done by imaginary time evolution. In practice we iterate the action of the evolution 
operator exp (— rfri/gp) (where ifgp is given in (|TB|)) on an initial wavefunction, renormalizing 
the wavefunction and calculating fj, at each time step dr, until /z reaches a stationary state 
within the required accuracy. If dr is small enough, we can approximate the evolution of the 
wavefunction in each interval dr by splitting the exponential exp {—dr Hgp) into the exponential 
of kinetic energy and the exponential of potential energy: 

0(r + dr) = exp {-dr Hgp) 0(r) ~ exp {-dr Edn) exp {-dr E^ot) 0(r) (64) 

where: 

^cin = |- (65) 

2m 

Epot = Uiz) + Ng\(P{z,T)\^. (66) 
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Figure 1: Spatial density of non condensed atoms in presence of a condensate in a ID harmonic 
trap, for: = 10^, = 30huj, gir, = 0.01h(ij/mY^'^. The corresponding chemical potential 
is /i = lA.lhuj. Full line: direct diagonalization. Points with error bars: Monte Carlo simulation 
with 2000 realizations. The length scale used is the harmonic oscillator length a^o = 



imu. 



The advantage of the splitting in ( p4D is that the action of and E^ot on a vector can be 
easily calculated by going back and forth in the Fourier space and in the real space. This 
procedure is called the splitting method and is it a common technique to solve the steady-state 
Gross- Pit aevskii equation. Let us now consider the Brownian motion simulation: because of 
our choice (|50| ) and ([5T| ) we have to calculate the action of the exponential of (3C/2 over a 
vector. It is convenient to write: 



exp 



/5 



C \ = [exp ((ir£)]" with 



dr 



/3/2 



n 



(67) 



The idea underlying equation (^3) is to choose n large enough so that the splitting approx- 
imation sketched above can be applied. Before implementing the splitting however, we still 
have to perform some transformation. We notice that if acting on functions orthogonal to the 
condensate wavef unctions, exp {drC) can be expressed as: 



exp (drC) 



Q 






Q* 



exp (drC 



9PJ 



(6J 



12 



where Cgp is the operator obtained by hnearization of the Gross-Pitaevskii equation, which 
is the same as C given by ([T5|) without the projectors. We are then led to the action of 
exp (drCgp) on a vector and will do it approximately by using the splitting technique. In order 
to be consistent however we must take care that in our simulation the approximate action of 
£gp corresponds exactly to the linearization of the approximate Gross-Pitaevskii equation in 
imaginary time that we have used to find the condensate wavefunction. We then obtain: 



exp {drCgp) ~ exp {-drEan) exp [-dr^E^ot - yu)] exp 



Ngdr 



(69) 



where the matrix in the last factor has the nice property of giving zero if squared. 

A last point concerns the choice of the time step dt of the Brownian motion, and of the 
total simulation time tmax guaranteeing that the steady-state is reached. 

The time step dt should satisfy ttmaxC^^ ^ 1 where amax is the largest eigenvalue of the 
friction matrix Eq.(BO). To estimate Omax we consider the eigenmode (wmax, '^max) of C with the 



largest eigenvalue Cmax- This value is well above the chemical potential yU, so that Vraax can be 
neglected as compared to Wmax, in the so-called Hatree-Fock approximation. In this case this 
largest energy eigenmode is an approximate eigenvector of the friction matrix, and we get the 
condition 

"maxC^i ^ o?^sinh(^emax)//5 < 1- (70) 

In a one- dimensional harmonic trap one expects Cmax — ~ in the Hartree-Fock approxi- 
mation, where e^^-^ is the maximal energy level of the bare trapping potential. 

The condition to be satisfied by the duration tmax of the simulation is amin^max ^ 1 where 
amin is the smallest eigenvalue of the friction matrix. We assume here that the corresponding 
eigenmode of a has components on energy modes of C with energy much smaller than fc^T only. 
In this case a^i^ corresponds approximately to the minimum eigenvalue of rjC In the case of 
real condensate wavefunction, the determination of the minimum eigenvalue of orthogonal 
to the condensate wavefunction boils down to the calculation of the first excited state energy 



tmaxT^f / 101 V) »1. (71) 



of Hgp — /i ; following we then find the condition 

2m 

In the weakly interacting regime one finds tmax ^ T^^- In the opposite Thomas- Fermi regime one 
finds that tmax scales as /i S> fiu). In practice, in the case of figure |l| we had: dt = 1.15xl0~^//ia; 
and tmax = 56.9 /hu. 



6 The time-Dependent Case 

Let us consider the following situation: initially we have a Bose-Einstein condensate in a 
thermodynamically stable state at thermal equilibrium. Suddenly a parameter of the system is 
changed (e.g. the trapping potential or the scattering length) so that the system is no more in 
an equilibrium state and it undergoes a dynamical evolution during a given time, after which 
we perform a measurement. 
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In order to study this situation we need to recall some more results of the time-dependent 
number conserving Bogolubov theory of j3^. It is shown in this paper that the expectation 
value of a one-particle observable X = J2i=i is given by: 

(X) ^ (iV-(5iV))(0|X(l)|0) + (0|X(l)|0(2)) + (0(2)|X(l)|0) 

+ J drj dr'{f\X (1)1?) {A\f)k{r')) + o(--^j . (72) 



In this equation the first term contains the leading contribution to {X) scaling as N in the 
limit (P), where is solution of the Gross-Pitaevskii equation (|I3|) , and 

{SN) = [ df{A\f,t)A{f,t)) , (73) 



so that N — (SN) is the mean number of particles in the condensate. The other terms involve 
the first correction to condensate wavefunction 0(2) and the contribution of the non condensed 
atoms to the expectation value of the observable. The equation of motion of 0(2) given in |Q 
is: 



where C is given by Eq.(|T5D and 

R{f) = -gN\(j){f)f(j){f){l + j dsA\s)A{s)) 

+ 2gN(j){f){A\f)A{f)) + gN(j)*{f){A{f)A{f)) 

- gN [ ds\ 0(s ) |^( [At(s )0(s ) + A(s )0*(s )1 A(f )) . (75) 



The first term in Eq.(|75D corrects the overestimation of the number of condensed particles in 
calculating their mutual interaction {N —>■ N — {1 + (SN)) in the Gross-Pitaevskii equation. 
The terms in the second line describe the interaction of the condensed particles and the non 



condensed ones. We refer to |22| and |33| for more details. 

The implementation of our method for the time-dependent case then consists of three parts: 
(i) we generate many stochastic realizations of the initial state sampling the correct equilibrium 
distribution, (ii) We perform a deterministic evolution in real time of the stochastic initial 
states, (iii) We calculate the expectation value by averaging over the stochastic realizations. 

Let us first concentrate on part (i): which is the condensate wavefunction to the lowest 
approximation, is found by solving the Gross-Pitaevskii equation (^) in steady-state. Once 
we have 0, we can use the Monte Carlo simulation described in section ^ to generate many 
realizations of the stochastic functions A and A* representing the non condensed fraction. To 
calculate 0(2), we calculate the source term ( ffS] ) and solve (|7^) at steady-state with an imaginary 
time evolution: 

' ^ ^^'^'^ ^ - -^vcit ^ 0) ( tr^ ) - ( j^i'rT.^ : S V (76) 



dr V </'(2)W ; ' ' ' V ^u^) J V Q*it = m*it = 0) 
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We note in this respect that the last term in (|7S|) for a given reahzation of A factorizes as A(r ) 
times an r -independent integral over s. 

The real time evolution part (ii) is more straightforward since it is deterministic. The 
wavefunction evolves according to ([I3|). Each stochastic realization of A, A* evolves according 



to ([T^) , and the corrections 0(2) , 0*2) evolve according to (^4D . The averages in the source term 
R involving A and A* must be calculated at each time step. 



7 Conclusions and Perspectives 



We have put forward a Monte Carlo method to sample numerically the Gaussian density op- 
erator for the non condensed modes obtained in the Bogolubov approach, once it is expressed 
in term of the Wigner quasi-distribution function. This method allows to implement Bogol- 
ubov theory both in steady-state and in the time-dependent case without having to diagonalize 
big matrices, which makes it readily scalable to two or three spatial dimensions. We plan to 
use this method to study thermodynamics and finite temperature dynamics of Bose-Einstein 
condensates with vortices. 



In this paper we have sticked to the Bogolubov theory as presented in P3|. In this frame, 
for the time-dependent case, the Gross-Pitaevskii condensate wavefunction, the non condensed 
fraction and the first correction to the condensate wavefunction must be first obtained in the 
initial equilibrium state and then they must be evolved separately. 

A natural development of the method beyond the Bogolubov approximation consists in 
a Monte Carlo simulation which samples the Wigner distribution for the whole field ip in 
equilibrium at temperature T followed by the real time evolution of the field according to the 
classical field equation: 

ihdtij = hoi^ + g\ij\^ij . (77) 



The classical field equation ([77|), which is formally identical to the Gross-Pitaevskii equation, 
corresponds in fact to the approximate evolution in real time of the Wigner quasi distribution 
function (truncated Wigner) . 

The classical equations of motion for the field ([77|) has the disadvantage of considering 
only stimulated processes. Consider for example two colliding condensates in three dimensions, 
corresponding to an initial field of the form 

'(/'(r ) = A {exp{—ikz) + exp{ikz)) . (78) 



By evolving this field according to ([77|) we will populate only modes that are plane waves of 
momentum hnk, n integer (e.g. n = 2 corresponds to the stimulated process k + k ^ —k + 2k). 
In reality we expect to have binary collisions producing emerging atoms along any spatial 
direction. At the level of ip this would require a broken translational symmetry in the x-y 
transverse plane. Fortunately, the advantage of being in the Wigner point of view is that all 
the modes are initially "filled" by quantum noise. In this case all processes can be stimulated 
and all possible symmetries of ip are broken. 

We acknowledge useful discussions with Cora Shlyapnikov, Yuri Kagan and Anthony Leggett 
in the early stage of this work. Laboratoire Kastler Brossel is a unite de recherche de I'Ecole 
normale superieure et de I'Universite Pierre et Marie Curie, associee au CNRS. 
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